draft of February 1, 2008 



The Parker Instability under a Linear Gravity 8 

Jongsoo Kim 1,2 ' 5 , Seung Soo Hong 2,6 , and Dongsu Ryu 3, ' 

ABSTRACT 

o\ . 

A linear stability analysis has been done to a magnetized disk under a linear gravity. 
We have reduced the linearized perturbation equations to a second-order differential 
equation which resembles the Schrodinger equation with the potential of a harmonic 
oscillator. Depending on the signs of energy and potential terms, eigensolutions can be 
classified into "continuum" and "discrete" families. When magnetic field is ignored, 
the continuum family is identified as the convective mode, while the discrete family 
as acoustic-gravity waves. If the effective adiabatic index 7 is less than unity, the 
former develops into the convective instability. When a magnetic field is included, the 
continuum and discrete families further branch into several solutions with different 
characters. The continuum family is divided into two modes: one is the original 
. Parker mode, which is a slow MHD mode modulated by the gravity, and the other is 

a stable Alfven mode. The Parker modes can be either stable or unstable depending 
on 7. When 7 is smaller than a critical value 7 cr , the Parker mode becomes unstable. 
The discrete family is divided into three modes: a stable fast MHD mode modulated 
by the gravity, a stable slow MHD mode modulated by the gravity, and an unstable 
mode which is also attributed to a slow MHD mode. The unstable discrete mode does 
not always exist. Even though the unstable discrete mode exists, the Parker mode 
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dominates it if the Parker mode is unstable. However, if 7 > 7 cr , the discrete mode 
could be the only unstable one. When 7 is equal 7 cr , the minimum growth time of the 
unstable discrete mode is 1.3 x 10 s years with a corresponding length scale of 2.4 kpc. 
It is suggestive that the corrugatory features seen in the Galaxy and external galaxies 
are related to the unstable discrete mode. 
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1. INTRODUCTION 

Parker (1966) demonstrated that the equilibrium state of interstellar gas, magnetic field, and 
cosmic-rays under the vertical gravity could be unstable to long wavelength perturbations along 
the direction of the equilibrium magnetic field. Since then, many authors have studied the details 
of this instability. For instance, two-dimensional, nonlinear equilibrium states resulting from 
the Parker instability under a uniform gravity were derived by Mouschovias (1974). The effects 
of rotation (Shu 1974; Zweibel & Kulsrud 1975; Foglizzo & Tagger 1994), non-uniform gravity 
(Horiuchi et al. 1988; Matsumoto et al. 1988; Giz & Shu 1993), and skewed magnetic field (Hanawa 
et al. 1992) on the Parker instability were discussed. The Parker instability was considered to 
be a plausible mechanism for the formation of large-scale interstellar clouds (Appenzeller 1974; 
Mouschovias et al. 1974; Blitz & Shu 1980; Shibata & Matsumoto 1991; Gomez de Castro & 
Pudritz 1992), a possible process for dynamo in accretion disks (Tout & Pringle 1992), and the 
force driving emergent flux tubes which are closely related to the active phenomena of the Sun 
(Shibata et al. 1989a, 1989b; Kaisig et al. 1990; Nozawa et al. 1992). 

The original analysis of the instability by Parker (1966) was performed under a uniform 
gravity, although he commented on the effects of a non-uniform, linear gravity. Since the uniform 
gravity has a jump at the mid-plane, the physical quantities should have a cusp. Because of 
this, the perturbations in the upper hemisphere cannot communicate with those in the lower 
hemisphere. So Parker's analysis should have been restricted to the "symmetric" mode where 
the vertical velocity perturbation vanishes at the mid-plane. In order to remove this difficulty, a 
non-uniform gravity which changes smoothly at the mid-plane should be introduced. 

Giz & Shu (1993) considered the Parker instability under a gravity whose functional form is 
—g t&nh.(z/2H*), where g Q and H* are constants, and z is the vertical distance from the mid-plane. 
With the scale height of stars for H*, this gravity closely represents the Galactic gravitational 
field. In the limit \z\ —>■ oo it is reduced to a uniform gravity, while in \z\ — ► it is reduced to a 
linear gravity. With this gravity, Giz & Shu obtained a dispersion relation which matches that 
under the uniform gravity. They called it the "continuum" family. But, in addition, they found 
a dispersion relation which has no correspondence in the uniform gravity. They called the latter 
the "discrete" family. Before Giz & Shu, Horiuchi et al. (1988) considered the Parker instability 
under a point-mass-dominated gravity, which is continuous across the mid-plane. It is a good 
model for the gravitational field in accretion disks. They considered isothermal perturbations 
with the effective adiabatic index 7 = 1 and obtained dispersion relations for both the symmetric 
and antisymmetric modes. Their solutions correspond to the continuum family by the Giz & 
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Shu classification. They did not see the solutions of the discrete family, since their analysis was 
restricted to isothermal perturbations. 

In this paper, we consider a linear gravity whose functional form is —g'z, where cf is a 
constant. The linear gravity represents fairly well the z-dependence of the Galactic gravity up to 
500 pc (Oort 1965), as well as the "vertical" gravity of accretion disks. It is also easier to handle 
than those introduced by Giz & Shu (1993) and Horiuchi et al. (1988). Although Parker (1966) 
considered a linear gravity, all he did was to show that the equilibrium state of an interstellar 
gas-field system becomes also unstable under it. Here, we investigate the resulting eigensolutions 
in detail: we classify them into the continuum and discrete families, discuss their nature, and give 
criteria for the existence of the unstable continuum and discrete families. Although the existence 
of the discrete family was already discussed by Giz & Shu (1993), our work would serve to get a 
better insight to its nature and its possible roles in astrophysics. 

Several authors have reported that there exist corrugations in our Galaxy (Quiroga & 
Schlosser 1977; Lockman 1977; Spicker & Feitzinger 1986; Alfaro et al. 1992; Malhotra 1994) and 
external galaxies (Florido et al. 1991) as well. Nelson (1976) suggested that the corrugations 
may be a manifestation of hydrodynamic acoustic-gravity waves. The difficulties with Nelson's 
suggestion are that the characteristic scale does not come out naturally and that the structures 
formed by the acoustic-gravity waves should be oscillatory and transient. Instead, in this paper 
we suggest that the corrugations may be formed as the results of the unstable discrete mode. 

The plan of the paper is as follows. In §2 the classification of the eigensolutions into the 
continuum and discrete families is given. In §3 the dispersion relations of the eigensolutions are 
analyzed in details. We first consider the dispersion relations of non-magnetized disks and then 
move to those of magnetized disks. In §4 a summary and discussion are given. The derivation of 
eigenvalue equations and the method we used to obtain the dispersion relation of the continuum 
family are described in appendices. 

2. CONTINUUM AND DISCRETE FAMILIES 

We consider the stability of a disk composed of gas, magnetic field, and cosmic-rays under 
a linear gravity. For the unperturbed equilibrium state, we follow Parker (1966) to assume that 
the gas is isothermal and that the ratio of magnetic pressure to gas pressure, a, and the ratio of 
cosmic-ray pressure to gas pressure, (3, are spatially constant. Then, the unperturbed gas density, 
gas pressure, cosmic-ray pressure, and magnetic pressure are given by a gaussian distribution in z. 
The equations for the stability analysis are obtained by introducing an infinitesimal perturbation 
to the unperturbed equilibrium state. Those equations can be combined into a single eigenvalue 
equation, which is a second-order differential equation. The detailed manipulation of the equations 
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is described in Appendix A. The resulting eigenvalue equation is 



where E and V Q are defined by 



+ {E- V ( 2 )i> = 0, (1) 



[l + a + p) + (2a + 7 )(^ 2 + u 2 ) ft 2 + 2a 7 (^ 2 + rfty 



Va 



2 (2a + 7) ft 2 - 2a 7 f 2 

1 2a(l + a + /3)i/ 2 ^ 2 (1 + a + + a + /? - 7 )(z/ 2 + z^ 2 ) 



4 (ft 2 - 2a^ 2 ) [(2a + 7 )ft 2 - 2a 7 ^ 2 ] (2a + 7 )ft 2 - 2a 7 i/ 2 ' 

Here £ is a normalized z-coordinate, v x and are normalized radial and azimuthal wavenumbers, 
$7 is a normalized angular frequency, and tp is a quantity proportional to the normalized vertical 



velocity multiplied by the square root of the unperturbed normalized gas density (see Eqs. [A.28] 



and [A30]). It should be pointed out that t/j 2 is the kinetic energy density associated with the 



motion in the z-direction. Normalization units are the disk scale height H, sound speed of 
unperturbed state a s , and unperturbed gas density at the mid-plane /0 o (0). Due to the simple 
functional form of the gravity we used, we could get an eigenvalue equation which resembles the 
Schrodinger equation with the potential of a harmonic oscillator. As the boundary condition, we 
assume that ip stays finite as £ — > 00. That is, the energy density perturbation does not diverge at 
infinite. 

The solutions of the eigenvalue equation can be classified into four different types by the signs 
of E and V Q (Figure 1). If E > and V Q < so that E is greater than V(() in the whole range 
of C, — 00 < ( < +00, all the positive values are allowed as the eigenvalue E (Type I). In this 
case, the eigenfunction if) oscillates as |£| increases. If E < and V D < so that E is smaller than 
V(C) in a finite interval which includes £ = and greater than V(() in the rest, all the negative 
values are allowed as the eigenvalue (Type II). The eigenfunction declines exponentially in the 
finite interval but stays oscillatory outside it. If E > and V Q > so that E is greater than V(() 
in a finite interval but smaller than V(C) m the rest, only discrete positive values are allowed as 
the eigenvalue (Type III). The eigenfunction oscillates in the interval but declines exponentially 
outside it. If E < and V Q > so that E is smaller than V(C) in the whole interval, there are no 
physically meaningful solutions satisfying the boundary condition (Type IV). With the tangent 
hyperbolic gravity, Giz & Shu (1993) also obtained an equation similar to Eq. (|]) and classified 
its solutions into the continuum and discrete families. According to the Giz & Shu's classification, 
the solutions of Type I and II correspond to the continuum family and those of Type III to the 
discrete family. 



2.1. Continuum Family 



For the continuum family with negative V Q , we introduce a new variable £ defined by 
£ = \/2 \/—V ( and transform Eq. (|l]) into that for the parabolic cylinder function (Abramowitz 
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& Stegun 1970): 



d£ 2 \4 



a = (3) 

As |f | increases, the parabolic cylinder function oscillates in the interval between nodal points. 
For |f | -> do, i) oc exp(±if 2 )/ 1 /f- So the energy density perturbation of the continuum family 
decreases as l/|f| at large |f|. It is interesting to notice that, for |f| — > oo, the z wavelength 
decreases as l/|f|, hence the group velocity along the z-direction increases as |f|. Consequently 
the energy flux of the continuum family stays constant at large |f |. 

The dispersion relations for the symmetric mode can be calculated by considering the lower 
boundary condition at the mid-plane, ip(£ = 0) = 0, while the dispersion relations for the 
antisymmetric mode by taking dip(^ = 0)/c£f = (see Appendix B). With the condition if)y/% = 
finite at |f | — > oo, there exists, for a given pair of wavenumbers u 2 and v 2 , an interval of Q 2 values. 
Here, we find solutions with the energy density perturbation zero at the upper boundary which 
is identified to the first nodal point from the mid-plane, f no de- Then, the value of Q 2 is uniquely 
determined for the given nodal point. Since the nodal point indicates the position where the 
vertical velocity goes to zero, the distance to this nodal point may be considered as an "effective" 
wavelength along the z-direction. A method to find the value of Q 2 with a chosen nodal point is 
described in Appendix B. 



2.2. Discrete Family 

For the discrete family with positive V Q and positive E, Eq. (|]) has exactly the same form as 
the Schrodinger equation with the potential of a harmonic oscillator. The solution of the harmonic 
oscillator problem is well known: the eigenfunctions are 



MO = 




and the corresponding eigenvalues are 

= 2n + l, (5) 

where H n is the Hermite polynomial of degree n, and n takes values 0, 1, 2, • • •. In this case, ip 
decreases exponentially at large |£|, so the energy density perturbation of the discrete family does 
as well. 

The eigenfunctions for n = 1, 3, • • • which have tp = at £ = represent the symmetric 
mode, while the eigenfunctions for n = 0, 2, • • • which have dip/dC, = at £ = represent the 
antisymmetric mode. The eigenvalue relation of Eq. (||) gives the dispersion relation, i.e., Q 2 as a 
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function of u 2 , Vy and n. Here, n indicates the number of nodal points along the z-direction, so it 
may be considered as the "effective" wavenumber in the z-direction. 



3. WAVES AND INSTABILITIES 

In this section, we explore the nature of the continuum and discrete families in detail by 
analyzing both the stable (waves) modes and the unstable (instability) modes. We first consider 
the case of non-magnetized disks, then include a magnetic field and cosmic-rays. 



3.1. Non-Magnetized Disks 

Without a magnetic field, the cosmic-rays are ignored, so a = (3 = 0. Without a magnetic 
field, there is no preferred direction in the x—y plane, and we can set v x = without losing any 
generality. Then, E and V Q become 

+ (6) 
2 7 



y o = 7 + n : y - (7) 



i (i-tK 

4 7^ 2 

Eq. (|l]) with the above E and V Q is exactly the one Nelson (1976) derived to study hydrodynamic 
acoustic-gravity waves under a linear gravity. 

In the case of 7 = 1 for isothermal perturbations, V a = 1/4 is always positive. So, Eq. (|l|) 
has only two types of solutions, Type III and Type IV, depending on the sign of E. The locus 
of E = is drawn with a dashed line in the phase space of v 2 and Q 2 in Figure 2a. We are 
interested only in the upper region of the line with E > where the solutions of Type III exist. 
The eigenfunctions of the solutions in this region are 

VKC) = ~^^=H n (4=) exp (-£) , (8) 



V2tt \/2 n n! VV2 
and the corresponding dispersion relations are 

n 2 = i/J + n+l. (9) 

Eq. (^) reduces to the dispersion relation f2 2 ~ v 2 . of sound waves for v 2 S> n + 1. So the dispersion 
relations are of the acoustic-gravity waves which have the characteristics of sound waves modulated 
by the vertical gravity, as was pointed out by Nelson (1976). In Figure 2b, the dispersion relations 
of the acoustic-gravity waves with n = 0, 1, • • • ,6 are plotted with solid lines. In this case, no 
convective mode exits, so all the modes are stable. 
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If 7 7^ 1) the signs of both E and VJ, can be positive or negative. In Figures 3a and 4a, three 
loci of E = 0, V = 0, and n 2 = are drawn for 7 = 0.8 and 1.3 with a dashed line {E = 0), a 
dotted line (V = 0), and a long-dashed line ($7 2 = 0). The phase space of v 2 and Jl 2 is divided into 
four regions by the three loci, which are identified as those of the Type II, III, and IV solutions. 
We are interested only in the regions of the Type II (continuum family) and the Type III (discrete 
family) solutions. In Figures 3b and 4b, we plot the dispersion relations of the discrete family 
solutions with n = 0, 1, • • • , 6 which are obtained from Eq. (||), and the dispersion relations of 
the continuum family solutions with nodal points at Cnode = 3, 5, 10 which are computed by the 
method described in Appendix B. 

The solutions of the discrete family, which are identified as acoustic-gravity waves, always 
have positive Q 2 , so they are stable. On the other hand, the solutions of the continuum family, 
which can be identified as the convective mode, are stable if 7 > 1, but unstable if 7 < 1, 
developing into the convective instability (Ryu & Goodman 1992). The mode with a larger Cnode 
is more unstable, because the gravity is a linearly increasing function of C- 

In the case of a uniform gravity, the square of the maximum growth rate of the convective 
instability is Q 2 = —(1 — 7^7 (Ryu & Goodman 1992). For 7 = 0.8, it is —0.25. Except at small 
wavenumbers, the growth rate of the convective instability under the linear gravity in Figure 3b is 
larger than the maximum growth rate under the uniform gravity. 

The dispersion relations of the continuum family solutions shown in Figures 3b and 4b are 
for the symmetric mode. In order to see the effects of the boundary condition at the mid-plane, 
we have also calculated the dispersion relations of the antisymmetric mode. Figure 5 shows the 
resulting dispersion relations with nodal points at Cnode = 2 and 3. Solid lines represent the 
dispersion relations of the symmetric mode, and dotted lines those of the antisymmetric mode. 
The growth rates of the both modes are nearly the same, unless Cnode 1- This is true, even if 
the magnetic field is included. The reason for this is that the gravity is small near the mid-plane. 
So the details around the mid-plane are not important unless Cnode 1- I n the rest of the paper, 
we then present only the dispersion relation of the symmetric mode for the continuum family 
solutions. 



3.2. Magnetized Disks 

3.2.1. Limit of H — ► 00 

Taking the limit H — > 00 means also making g' — > 0. So the problem reduces to that of finding 
the dispersion relation in a uniform medium without gravity. Since the perturbation of cosmic-ray 
pressure is zero in the uniform medium (see | A26| ), we should get the dispersion relations for three 
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magnetohydrodynamic (MHD) waves. From ( |A18|) to ( |A25D , by setting d/dz = 0, we have 



2 ^cxky 



UJ 



— - (2a + 7 )(A£ + k 2 y ) — + 2a 7 (^ + 



^ = 0, 



(10) 



where a s is the isothermal sound speed. The above dispersion relation can be decoupled into one 
for the Alfven waves 

(11) 



u? - v\k 2 = 0, 



u 

and the other for the slow and the fast MHD waves 

oj 4 - (v\ + c 2 s )(k 2 x + k 2 y )u 2 + v 2 A c 2 k 2 y (k 2 x + k 2 y ) = 0. 
Here, va = 2aa 2 and c 2 = 7a 2 are the Alfven and adiabatic sound speeds. 



(12) 



3.2.2. Case of i/ x = 

It is known that the growth rate of the Parker instability is larger for perturbations with 
larger x wavenumbers (Parker 1967). However, in this subsection we first consider the case of 
v x = in detail, since the equations involved are much simple and yet most important physical 
aspects are included. We comment briefly on the general case of v x ^ in the next subsection. 

In the case of v x = 0, the x-components of velocity and magnetic field perturbations are 
decoupled from the rest and result in Alfven waves with the dispersion relation (pd|) (see Eqs. [ A19| 



and [A22|). So, with the rest of the perturbation equations, we get the eigenvalue equation (|l]) 



with E and V Q given by 



(2a + 7)^ 2 -2«7i/2 



1 (l + a + /?)(l + a + /?-7)^ , . 

° 4 (2a + 7)ft 2 - 2a 7 ^2 ' 1 ' 

In order to classify its solutions, we draw in Figure 6a the loci of E = 0, V a = 0, and 
(2a + 7)J7 2 — 2a^v 2 = with two dashed lines, a dotted line, and a long-dashed line, respectively, 
for a = 1, (3 = 0, and 7 = 1. Since E = reduces to a quadratic equation of O 2 , there exist 
two branches, which are labeled by E + = (the positive branch) and E~ = (the negative 
branch). The loci of E~ = and (2a + 7)£1 2 — 2a7f 2 = intersect at the point where 
v 2 = (2a + 7)(1 + a + f3)/(2aj) and J7 2 = 1 + a + [3. The intersection point exists also in the 
case of v 2 7^ 0, and its position is determined only by a, (5, and 7. So, its location for the case 
of v 2 7^ is the same as that for the case of v 2 = 0. The four loci divide the phase space of v 2 
and f2 2 into six regions, which are identified in Figure 6a by the type of solutions. The continuum 
family solutions exist in the regions of Type I and II, and the discrete family solutions in Type 
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III. Figure 6b shows the dispersion relations of the discrete family solutions with n = 0, 1, 2 and 
those of the continuum family solutions with nodal points at Cnode = 3, 5, 10. 

The continuum family solutions correspond to those of the original Parker instability and are 
identified as a slow MHD mode modulated by the gravity (Parker 1966; Shu 1974). For the same 
nodal point and the same parameters (a, /3, 7), the growth rate of the continuum family is larger 
under a linear gravity than under a uniform gravity. A full comparative study of the growth rate 
under different models of gravity will be given in a separate paper. 

The discrete family solutions in the upper left region of Figure 6b are MHD counterparts of 
Nelson's acoustic-gravity waves, which are identified as fast MHD waves modulated by gravity. The 
discrete family solutions in the region enclosed by the loci of E~ = and (2a + 7)f2 2 — 2a^v 2 = 
represent a new mode which was not discovered in any previous work. They are always stable, 
since Q 2 > at the intersection point of the loci. On the following grounds, we identify them as 
slow MHD waves modulated by gravity. For v y — > 00 and n — > 00, the dispersion relation of the 
discrete family becomes 

ft 4 - (2a + 7 )K 2 + vl)® 2 + 2a 7 (^ 2 + v 2 z )v 2 y = 0, (15) 



if we set 2y/Vc,n = v z (see Eqs. ||, and |i4}] ). Note that V a becomes constant in this limit. 



Eq. (15) is the dispersion relation for the fast and slow waves (see Eq. [12]). The fast waves 



correspond to the solutions in the upper left region, while the slow waves to those in the region 
enclosed by loci of E~ = and (2a + 7)f2 2 — 2a^v 2 = 0. 

If the slope of the locus of V Q = in the phase space of v 2 and fi 2 is equal to or greater than 
zero, the continuum family solutions are stabilized. For given a and j3, the critical value of the 
adiabatic index, 7 cr , which gives zero slope, is 



l + %a + P 



7cr = — 3 — ■ ( 16 ) 



Note that this is same as the critical value for the case of a uniform gravity (Parker 1966). To see 
the nature of the solutions in the case of 7 > 7 cr , we draw in Figure 7a the loci of E = 0, V = 0, 
and (2a + 7)^ 2 — 2a^v 2 = for a = 1, (3 = 0, and 7 = 7 cr . Figure 7b presents the dispersion 
relations of the discrete family solutions with n = 0, 1, 2 and those of the continuum family 
solutions with nodal points at Cnode = 3, 5, 10. As in the case of Figure 6, there are solutions 
of the continuum family in the regions of Type I and II, but now they are stable with Q 2 > 0. 
Also there are solutions of two stable discrete modes in the upper left region and in the region 
enclosed by loci of E~ = and (2a + 7)f2 2 — 2a7^ 2 = 0. However, there is another region of 
Type III enclosed by loci of E~ = and V = 0, as can be seen in the enlarged plot of Figure 7c. 
The solutions of the discrete mode in this region are unstable with Q 2 < 0. There are an infinite 
number of unstable solutions with n = 0, 1, 2, • • -, and Figure 7d shows the dispersion relations 
with n = 0, 1, 2. 
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Figure 8 shows the eigenfunctions of the unstable discrete mode solutions with n = 0, 1, 2 for 
a = 1, f3 = 0, 7 = 7 cr , and v x = 0. For a given n, the most unstable growth rate and its azimuthal 
wavenumber are used. The eigenfunctions of odd n's are of the symmetry mode with ip(Q) = 0, 
while those of even n's the antisymmetry mode with dip(0)/d( = 0. As was noted in §2.2, the 
energy density perturbation of the discrete mode solutions is concentrated around C = 0. Since 
the stable discrete mode solutions in the upper left region are identified as the fast MHD mode 
modulated by gravity, we may attribute the unstable discrete solutions to the slow MHD mode 
modulated by gravity as the stable discrete solutions in the region enclosed by loci of E~ = and 
(2a + 7 )f] 2 - 2a-/u 2 = 0. 

The unstable discrete mode does not always exist. Its existence depends on whether the two 
loci of E~ = and V Q = meet at two points: one at the origin {y 2 = 0, O, 2 = 0) and the other 
denoted by (z^ )int , ^ 2 nt )- For f3 = 0, Figure 9a shows the equi-f 2 int contours on the [a, 7] plane, 
and Figure 9b does the equi-f2 2 nt contours. The unstable discrete mode exists if 7 has a value 
in the region between the two lines of flf nt = in Figure 9b, or if 1 < 7 < 7i nt . Here, 7i nt is a 
function of a in this case of = and u x = 0, but it depends on other parameters, too, in general. 
If a medium has 1 < 7 < min(7 cr , 7i n t), both the unstable continuum and discrete modes exist. 
The growth rate of the unstable continuum mode, especially at small wavelengths, is always larger 
than that of the unstable discrete mode, making the continuum mode dominate over the discrete 
mode. On the other hand, if 7 cr < 7 < 7i nt , the discrete mode is the only unstable one, so it 
determines the evolution of the medium. 

If the scale height (160 pc) and the velocity dispersion (6.4 km/sec) of interstellar clouds are 
substituted for H and a s , respectively, the growth time of the most unstable discrete solution with 
n = becomes 1.3 x 10 8 years at the wavelength of 2.4 kpc. 

3.2.3. Caseofv x ^0 

When v x / 0, the Alfven mode is no longer decoupled from the other modes. Figure 10a 
presents the loci of E = 0, V a = 0, (2a + j)n 2 - 2ajv 2 = 0, and ft 2 - 2av 2 = for a = 1, (3 = 0, 
7 = 1, and v x = 1. In this general case, both loci of E = and V Q = are divided into positive and 
negative branches. The locus of — 2av 2 = is drawn with a dot-dashed line. The six loci divide 
the phase space of v 2 and Q 2 into eight regions. Since 7 = 1 there is no region containing the 
unstable discrete mode solutions. The solutions of the continuum and discrete families identified 
in Figures 6 and 7 are also found in the figure, and Figure 10b shows their dispersion relations. 
In addition, we have a region of Type II enclosed by the loci of = and £l 2 — 2av 2 = 0. The 
solutions in this region are of the continuum family and always stable. Since one of the boundaries 
of this region corresponds to the dispersion relation of the Alfven waves, we attribute its solutions 
to the Alfven mode. 
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4. SUMMARY AND DISCUSSION 

We have performed a linear stability analysis for magnetized disks under a linear gravity. The 
linearized perturbation equations have been reduced to an eigenvalue equation that has the same 
form as the Schrodinger equation with an energy E and a potential V D C 2 , where E and V Q are give 
by Eqs. (|A3^) and (|A33|) . By the signs of E and V , the eigensolutions have been classified into 



two families: the continuum family for V a < (independent of the sign of E) and the discrete 
family for E > and V Q > 0. 

Without a magnetic field, the solutions of the continuum family are identified as the convective 
mode and those of the discrete family as the acoustic-gravity waves which were previously discussed 
by Nelson (1976). The convective mode develops into an instability if 7 < 1 (Ryu & Goodman 
1992). 

With a magnetic field, the solutions of the continuum family are divided into two modes: 
a stable Alfven mode, and the original Parker mode. The Parker mode is a slow MHD mode 
modulated by the gravity, and becomes unstable when the effective adiabatic index is smaller 
than the critical value 7 cr = (1 + a + (5) 2 / '(1 + 3a/ 2 + (3) for v x = 0. The solutions of the discrete 
family are divided into three modes: a stable fast MHD mode modulated by the gravity, a stable 
slow MHD mode modulated by the gravity, and an unstable mode. The two stable discrete modes 
always exist, but the unstable discrete mode exists only if 1 < 7 < 7; n t, where 7; n t is a function 
of the parameters of the problem. The unstable discrete mode was first discovered by Giz &; Shu 
(1993). 

In the Galactic context, we discuss the possible roles of the unstable modes of the continuum 
and discrete solutions on the formation of giant molecular clouds and corrugations. For the 
scale height and the velocity dispersion of the cloud distribution, we take 160 pc and 6.4 
km/sec (Falgarone &: Lequeux 1973). Then, the maximum growth time and the corresponding 
length of the continuum solution are 1.5 x 10 7 years and 340 pc, for the parameters of 
a = 0.25,/? = 0.4,7 = 1.0, v x = 0.0 and Cnode = 3. Here, the values of a and /3 are the canonical 
ones of the interstellar medium (Spitzer 1978) and the chosen nodal point is ~ 500 pc (Elmegreen 
1982). On the base of the time and length scales of the Parker instability under the uniform 
gravity, Mouschovias et al. (1974) preferred the interstellar medium behind Galactic shocks to the 
canonical interstellar medium since the growth time of the Parker instability with the canonical 
parameters of the interstellar medium is larger than the life time of the giant molecular clouds, 
3 x 10 7 years (Blitz & Shu 1980). Under the linear gravity, the unstable modes of the continuum 
solution may, however, develop within the life time. Therefore, developments of the unstable 
modes of the continuum solution of the Parker instability under the linear gravity are not confined 
to spiral arms but ubiquitous in the Galaxy. 

We also suggest that the corrugations observed in our Galaxy and external galaxies may 
have been formed by the unstable discrete mode. For parameters, a = 1, /3 = 0, 7 = 7 cr , and 
v x = 0, the growth time of the most unstable discrete solution with n = is 1.3 x 10 8 years at the 
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wavelength of 2.4 kpc. This time scale is comfortably shorter than the age of the galaxies, and 
the length scale is consistent with the characteristic wavelength of the vertical structure of the 
Carina-Sagittarius spiral arm (Alfaro et al. 1992). Furthermore, the antisymmetric nature of the 
n = solution is also consistent with the corrugations found by Alfaro et al. (1992). 

In this paper, we have neglected the effects of rotation on the Parker instability. Although the 
effects of differential rotation which occurs in accretion disks as well as galactic disks are difficult 
to work with analytically, those of rigid-body rotation can be included in the frame of the present 
work. The inclusion of rigid-body rotation is expected to have two major effects on our results. 
First, although it does not change the criterion of the instability, the inclusion of rigid-body 
rotation reduces the growth rate as already pointed in Shu (1974) and Zweibel & Kulsrud (1975). 
Second, also as pointed in Shu (1974), with rigid-body rotation the Lindblad oscillation enters to 
the problem making the classification, for instance the one in Figure 10, more complicated. We 
are working on these and the results will be reported in a separate paper. 

At an early stage of this work JK was supported by DAEWOO Post-Graduate Scholarship. 
The work by SSH and JK was supported by Seoul National University DAEWOO Research Fund. 
The work by DR was supported in part by Seoam Scholarship Foundation. We are grateful to 
Dr. T. W. Jones and the referee, Dr. E. G. Zweibel, for comments on the manuscript. 



A. EQUILIBRIUM STATES AND LINEARIZED PERTURBATION EQUATIONS 

The MHD equations for gas with cosmic-rays are given by 

^+V-(pv) = 0, (Al) 

— + v ■ VvJ = - Vp - Vp cr - V— + —B-VB + pg, (A2) 
^=Vx(vxB), (A3) 

B ■ V Pcr = 0, (A5) 

where p cr is the cosmic ray pressure. 

For cosmic-rays we use the equation which was employed in the study of the effect of rotation 
on the Parker instability by Shu (1974). It is valid on the assumption that cosmic-rays adjust 
themselves to flows with an infinite speed along magnetic field lines. If the detail dynamics of 
cosmic-rays is important, for instance, around the shocks where they are accelerated, it can be 
treated either by the diffusion-advection equation such as derived by Skilling (1975) or by the 
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two- fluid model (Drury & Volk 1981). In the two-fluid model, cosmic-rays are represented by 
a fluid with an adiabatic index between 4/3 and 5/3. With these more realistic treatments 
of cosmic-rays, it was shown that fluids with cosmic-rays are subject to cosmic-ray induced 
instabilities with a scale comparable to the diffusion length of cosmic-rays (see, e.g. Drury & Falle 
1986; Zank, Axford, & McKenzie 1990; Kang, Jones, & Ryu 1992; Ryu, Kang, & Jones 1993; 
Begelman & Zweibel 1994). However, since here we are interested in the Parker instability which 
occurs in a much larger scale mainly by the coupling of magnetic field with gravity, the above 
simpler equation for cosmic-rays should be enough. 

To describe the local behavior of the Parker instability in magnetized disks, we introduce 
Cartesian coordinates (x, y, z) whose directions are radial, azimuthal and vertical, respectively. 
We assume that the gravity, g, is vertical and a linear function of z, g = (0, 0, — g'z), and 
that unperturbed magnetic field, B Q , is a function of z and has only an azimuthal component, 
B a = [0, B (z), 0]. Here, g' is a positive constant, We further assume that unperturbed gas is 
isothermal and the ratio of its magnetic pressure to gas pressure, a, and the ratio of its cosmic-ray 
pressure to gas pressure, /3, are spatially constant, i.e., p Q = a 2 p Q , B 2 /(8ir) = ap D , and p CTt0 = (3p , 
where a s is the isothermal sound speed, the quantities with subscript indicate unperturbed 
quantities. Then, the vertical distribution of unperturbed quantities is 



PoO) = Po(z) = Pcr,o( Z ) = B o( Z ) 

Po(0) p o (0) p cr , o (0) £2(0) 



exp 



1(± 

2 \H 



(A6) 



where H 2 = (1 + a + (3)a 2 s /g' and p o (0), p o (0), £> C r,o(0), B o (0) are the values of p (z), p (z), Pct,o( z ) 
B (z) at z = 0. 

We perturb the above equilibrium state. The perturbed state is described by 



v; p = Po + Sp; p = p + 5p; p cr = p CTfi + 5p CI ; B = B Q e y + SB. 



(A7) 



If we insert the quantities in Eq. ( |A7| ) into Eqs. ( |Al| )-(|Aq) and keep only the first-order terms of 
the perturbed variables, we obtain the linearized perturbation equations 



ddp p 



ZV Z + p 



dt H 2 ' 
dv x d5p d5p c 



/5 ° dt dx dx 



+ 



dv x dvy_ Qvz 
dx dy dz 

B d5B y B 86B X 

4tt dx 4tt dy 



0, 



0. 



dv y 85p d5p c 



^° dt dy dy 



+ 



B n 



8ttH 2 



z5B z = 0, 



dv z d6p d5p c - 
^° dt dz dz 



B Q B d5B z B ddBy (l + a + P)a 2 



8ttH 2 ' 



y 4vr dy 
d8B x 



4tt dz 



H 2 



z5p = 0, 



dt 



dy 



(A8) 

(A9) 
(A10) 
(All) 
(A12) 
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85 B v „ dv 



Bo 



dv z 



°8x 2H2 ZVz + B °dz 



d5B z dv z 



01 



0. 



dSp _ po_ 

at h 2 



d5p c 



dy 

dv x dvy_ dv : 
dx dy dz 

Pcr,o 



-5T ~ ~7F? ZVz + 7Po ( ^f 1 + + ) =0, 



0. 



(A13) 
(A14) 

(A15) 

at * — (A16) 

Since the coefficients of Eqs. (A8)-( A16|) do not depend explicitly on x, y, t, the perturbed 
quantities can be Fourier-decomposed with respect to x, y, t. So the perturbed quantities have 
the following form 

SQ(x, y, z; t) = SQ(z) exp(iuit — ik x x — ik y y), (A17) 

where u is the angular frequency and k x and k y are the wavenumbers along the x and y directions. 
Inserting the above into Eqs. (A8)-( A16 ), we get 



Po 



dv ? 



iojSp - Jp zv z + Po [ -ik x v x - ikyVy + — 



dz 



0. 



lUPoVx 



Bq Bq 

ik x 5p - ik x 5p cr - ik x —5B y + ik y —5B x = 0, 



iujp Vy - ikybp - iky5p CT + g7r ^2 zSBz = °> 



dbp ddp cl 
iujp v z + — — + — — 
dz dz 



B n 



8nH 2 



zbBy + iky—^-8B z + 

47T 



B d5B y (1 + a + (3)a 2 



47r dz 



+ 



# 2 



icu5B x + ik y B Q v x = 0, 



iuSB v — ik x B Q v x 



B n 



dv z 



2Hi zv z + B„- 



0, 



iu>5B z + ik y B v z = 0, 



tf 2 ' 



iwop - -775- zv z + 7p D -ik x v x - ik y v y + — — 



flz 



icaSp, 



Pcr,o 
cr rr 9 Zl> 2 



# 2 



(A18) 

(A19) 
(A20) 

z5p = 0, (A21) 
(A22) 

(A23) 
(A24) 

(A25) 
(A26) 



The above nine equations (Eqs. |A18| - flA26|l ) with nine unknowns (Sp, v x , v y , v z , 5B X , 5B y , 
5B Z , Sp, 5p cr ) can be combined into a single second-order differential equation for one unknown, 



d?v z 



dz 2 H 2 dz 



+ 



z dVz (l + a + p)/H 2 + (2a + 1 )(k 2 + k 2 ) to 2 /a 2 + 2a 7 (/c 2 + k 2 )k 2 

H — ; o , o = — V z 



(2a + 7 )w 2 /a 2 -2a-fk 2 



2a(l + a + l3)kluo 2 /a 2 s 



("7 



at — 2afc 2 



(2a + 7 )u; 2 /a 2 - 2a 7 fc 2 



(1 + a + + a + [5 - 7 )(/c 2 + k 
(2a + j)uj 2 /a 2 - 2a 7 /c 2 



z 2 

x w v * 



o. 



(A27) 
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With \& defined as 



v z = ^ exp 



1 / z x 2 
4 V# 



(A28) 



the equation reduces to a second-order differential equation without the first-order derivative term 



dPm J 1 "Vaf - [(1 + a + + (2a + 7 )(fc 2 + gj u^/af + 2a 7 (fc 2 + 

2^2 + (2a + -f)cu 2 /a 2 - 2a 7 fc 2 * 

2a(l + a + /3)£; 2 u;7a 2 (1 + a + + a + (3 - 7 )(A£ + k 



4 (w 2 /a 2 - 2afc 2 ) [(2a + 7 )w 2 /a 2 - 2a 7 fc 2 ] ( 2 « + 7)^ 2 /«s - 2a 7 /c 2 



x^ = 0. 



We define the dimensionless variables by 

Q = u;H/a s , v x = k x H, u y = k y H, ijj = *5?/a s , C = z/H. 
Then, the combined perturbation equation becomes 

d 2 v> 



+ {E- VoC 1 )^ = 0, 



(A29) 



(A30) 



(A31) 



where E and V Q are given by 



I Q 4 
E ^2 + - 



[l + a + P) + (2a + 7 )K 2 + v 2 ) n 2 + 2a 7 (z, 2 + u 2 )u 2 



(2a + 7 )fi 2 -2a 7 f 2 



2a(l + a + f3)u 2 Q 2 



4 (ft 2 -2az^ 2 ) (2a + 7 )0 2 - 2a 7 z/ 2 



+ 



(1 + a + /3)(1 + a + P - j)(u 2 + v 2 y ) 
(2a + 7 )ft 2 - 2a 7 f 2 



(A32) 
(A33) 



B. METHOD TO COMPUTE THE DISPERSION RELATIONS OF THE 

CONTINUUM FAMILY SOLUTIONS 



For given a, (3, and 7 , E and V Q in Eqs. ( A32| ) and ( |A33 ) are functions of v%,,Vy, and f2 2 . So 
getting the dispersion relations is equivalent to finding the values of v 2 , v 2 , and SI 2 which satisfy 



the boundary conditions imposed on Eq. ( |A31 ). We consider the upper boundary condition ip = 



at ( = Cnode and the lower boundary condition tp(0) = or dip(0)/dC = at the mid-plane. Here, 
Cnode is the nearest nodal point from the mid-plane. The condition ip(0) = doesn't allow gas 
to cross the mid-plane. We call the solutions that satisfy this condition as the symmetric mode. 
On the other hand, the condition difj(0)/d^ = allows gas to cross the mid-plane. We call the 
solutions that satisfy this second condition as the antisymmetric mode. 
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We introduce an index, i, and define 



(i = (i-l)A(, AC 



Cnode 

N-l' 



(Bl) 



where N is the total number of points in the interval, [0, Cnode]- With a central difference scheme 
the finite-difference representation of the perturbation equation ( |A3l| ) is 



^-i + [A( 2 (E - VoCf) - 2]Vi + = 0, 



(B2) 



where i runs from 3 to ./V — 2 for the symmetric mode and from 2 to N — 2 for the antisymmetric 
mode. Because ipN = at the upper boundary, the finite-difference equation at the point i = N — 1 
is 

+ [A( 2 {E - KC^-i) - 2]^tv-i = 0. (B3) 
For the symmetric mode with tpx = 0, we use the finite-difference equation at the point of i = 2, 

[A( 2 (E - VoCl) - 2]^ 2 + = 0. (B4) 

For the antisymmetric mode with i/jq = ip2, w e use the finite-difference equation at the point of 
i = l, 

[A( 2 (E-V ( 2 )-2}ij 1 + 2^2 = 0. (B5) 



We can write the simultaneous equations for ipi in a matrix form, 



Mip 



( * * 

* * 



\ 



\ ( i>s \ 



lpN-2 
V tpN-1 J 



( o \ 




V o / 



(B6) 



where s is 2 for the symmetric mode and 1 for the antisymmetric mode. Non-trivial solutions 
exist, only if the determinant of the tri-diagonal matrix M is equal to zero. For given a, (3, and 
7, the coefficients of the matrix M are functions of v 2 , v 2 , and So, the problem becomes an 
eigenvalue problem, i.e., with given v 2 and u 2 , we should find Q 2 which makes the determinant of 
M equal to zero. We used the Gauss-Seidel method for it. 
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Fig. 1. — By the signs of E and V a , the solutions of the eigenvalue equation (||) are classified into 
four types. For Type I and II, the solutions belong to the continuum family. For Type III the 
solutions belong to the discrete family. For Type IV, the solutions do not satisfy the imposed 
boundary condition. 

Fig. 2. — (a) Division of the phase space of v 2 and Q, 2 into regions with the solutions of different 
types for a non-magnetized disk with 7 = 1. The dashed line is the locus of E = 0. (b) Dispersion 
relations of the discrete family solutions with n = 0, 1, 2, • • • , 6 in the region of Type III. 

Fig. 3. — (a) Division of the phase space of v 2 and Q, 2 into regions with the solutions of different 
types for a non-magnetized disk with 7 = 0.8. Three loci of E = (the dashed line), V Q = (the 
dotted line), and S7 2 = (the long-dashed line) are drawn, (b) Dispersion relations of the discrete 
family solutions with n = 0, 1, 2, • • • , 6 in the region of Type III and the dispersion relations of 
the continuum family solutions with nodal points Cnode = 3, 5, 10 in the region of Type II. 

Fig. 4. — (a) Division of the phase space of v 2 and Q 2 into regions with the solutions of different 
types for a non-magnetized disk with 7 = 1.3. Three loci of E = (the dashed line), = (the 
dotted line), and fi 2 = (the long-dashed line) are drawn, (b) Dispersion relations of the discrete 
family solutions with n = 0, 1, 2, • ■ ■ , 6 in the region of Type III and the dispersion relations of 
the continuum family solutions with nodal points Cnode = 3, 5, 10 in the region of Type II. 

Fig. 5. — Dispersion relations of the symmetric mode solutions (solid lines) and the antisymmetric 
mode solutions (dotted lines) of the continuum family with nodal points Cnode = 2 and 3 for a 
non-magnetized disk with 7 = 0.8. 

Fig. 6. — (a) Division of the phase space of v 2 and Q 2 into regions with the solutions of different 
types for a magnetized disk with a = 1, /3 = 0, 7 = 1, and u x = 0. The loci of E = (the two dashed 
lines), V = (the dotted line), and (2a + 7) SI 2 — 2a^v 2 = (the long-dashed line) are drawn, (b) 
Dispersion relations of the discrete family solutions with n = 0, 1, 2 in the regions of Type III. 
The dispersion relations of the continuum family solutions with nodal points Cnode = 3, 5, 10 in 
the regions of Type I and II. 

Fig. 7. — (a) Division of the phase space of v 2 and Q, 2 into regions with the solutions of different 
types for a magnetized disk with a = 1, f3 = 0, 7 = 7 cr , and u x = 0. The loci of E = (the two 
dashed lines), V a = (the dotted line), and (2a+j)Q 2 — 2wyv 2 = (the long-dashed line) are drawn, 
(b) Dispersion relations of the discrete family solutions with n = 0, 1, 2 in the regions of Type III 
and the dispersion relations of the continuum family solutions with nodal points Cnode = 3, 5, 10 in 
the regions of Type I and II. (c) Magnification of the region of Type III enclosed by loci of E~ = 
and V = 0. (d) Dispersion relations of the discrete family solutions with n = 0, 1 and 2 in the 
magnified region of Type III. 

Fig. 8. — Eigenfunctions of the unstable discrete mode solutions with n = 0, 1, 2 for a = 1, /3 = 0, 
7 = 7cn and v x = 0. For a given n, the most unstable growth rate and its azimuthal wavenumber 
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are used. 

Fig. 9. — Contours of equi-f 2 int (a) and equi-f2 2 nt (b) for (3 = on the (a, 7) plane. The point 
(V^ int , fifnt) is defined by the intersection of the two loci of E~ = and V = 0. The intersection 
point exists only if 7 has a value in the region enclosed by the lines labeled with and 00. 

Fig. 10. — (a) Division of the phase space of f 2 and f2 2 into regions with the solutions of different 
types for a magnetized disk with a = 1, = 0, 7 = 1, and v 2 = 1. The loci of E = (the 
two dashed lines), V Q = (two dotted lines), (2a + j)Q 2 — 2a^Vy = (the long-dashed line), and 
Q? — 2aVy = (the dot-dashed line) are drawn, (b) Dispersion relations of the discrete family 
solutions with n = 0, 1, 2 in the regions of Type III and the dispersion relations of the continuum 
family solutions with nodal points Cnode = 3, 5, 10 in the regions of Type I and II. 
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